Mathematical Biosciences

Immune cell engagers (ICE) such as bispecific antibodies (bsAbs), within an immunological synapse, bind and link CD3 on a T cell to a target antigen (TAA) on a cancer cell, forming a trimer (CD3:bsAb:TAA complex). With sufficient trimer numbers within the synapse, the T cell can become activated and promote cancer cell killing. Elranatamab, a CD3-bispecific antibody for multiple myeloma, has received FDA and EMA filing acceptance (August 2023 and December 2023, respectively) adding to a growing list of bsAbs that are treating patients. In the drug development stages of ICE bsAbs, mechanistic modeling approaches are often used to attain a greater quantitative understanding of the modality, preclinically, and provide human pharmacokinetic and efficacious dose predictions to aide in Phase 1 trial design. To date, the majority of ordinary differential equation (ODE) trimer models treat the tumor compartment as well-mixed and trimer formation is governed by a bulk population reaction not accounting for individual synapses. This lack of discrimination can lead to imprecise analysis when analyzing results across E:T ratios using metrics like trimers per T cell or trimers per target cell. To this end we developed an ODE trimer model based on single-synapse complexes (one target cell/one immune cell) with 2D cross-linking trimer formation. We show computationally that the number of trimers per synapse is invariant to the value of the E:T ratio for a given free bsAb concentration, a property that cannot be captured by non-synapse models. A simple demonstration of this discrepancy using the well-known Betts trimer model is presented. We then apply the Betts trimer model coupled to a tumor growth inhibition (TGI) module to show that our synapse-based trimer model is easy to substitute in to model TGI, including the addition of a trimer-per-synapse activation threshold function for cell killing. Overall, our model attempts to balance mechanistic fidelity while limiting the complexity of the model.

There is a large variability between individuals in the response to anesthetic drugs, that seriously limits the achievable performance of closed-loop controlled drug dosing. Full individualization of patient models based on early clinical response data has been suggested as a means to improve performance with maintained robustness (safety). We use estimation theoretic analysis and realization theory to characterize practical identifiability of the standard pharmacological model structure from anesthetic induction phase data and conclude that such approaches are not practically feasible.

Determining optimal antibiotic dosing strategies is complex. Clinically, some antibiotics work best in continuous low doses, while others require high repeated pulses. However, a rational understanding of the best approach depending on the specific pairing of antibiotic and bacterial species remains unclear. Using mathematical models, we analyze bacterial populations under two strategies - constant concentration and repeated dosing - given fixed pharmacodynamic and pharmacokinetic properties. Our results reveal that the shape of the dose-response curve, which measures bacterial net growth rate against antibiotic concentration, is crucial. Specifically, its concavity determines the best strategy. In cases where the curve exhibits multiple concavities, additional factors such as tolerable dosing range influence the regimen. These findings challenge the universal application of "hit hard and hit early," as some recommended schedules include lower, constant doses. This work contributes to the literature on rational antibiotic prescription, aiming to minimize antibiotic use and combat antimicrobial resistance.

The interaction and possibly interference between viruses infecting a host population is addressed in this work. We model two viral diseases with a similar transmission mechanism and for which a vaccine exists. The vaccine is characterized by its coverage, induced temporary immunity, and efficacy. The population dynamics of both diseases consider infected individuals of each illness and hosts susceptible to one but recovered from the other. We do not incorporate co-infection. Two main transmission factors affecting the effective contact rates are postulated: i) the virus with a higher reproduction number can superinfect the one with a lower reproduction number, and ii) there exists some induced (indirect) protection induced by vaccination against the weaker virus that reduces the probability of infection by the stronger virus. Our results indicate that coexistence of the viruses is possible in the long term, even considering the absence of superinfection. Influenza and SARS-CoV-2 are employed to exemplify this last point, observing that the time-dependent effective contact rate may induce either alternating outbreaks of each disease or synchronous outbreaks. Finally, for a particular parameter range, a backward bifurcation has been observed for dynamics without vaccination.

Drug-target occupancy time--the cumulative duration a target remains bound--critically influences therapeutic efficacy. While Copelands widely-used residence time (1/koff) emphasizes dissociation kinetics, it neglects association rates, rebinding events, and drug elimination that affect in vivo outcomes. Returning to Paul Ehrlichs 1913 principle that drugs act only when bound ("Corpora non agunt nisi fixata"), we develop a mathematically rigorous framework defining Ehrlich occupancy time as the integral of fractional target occupancy over time. Our approach explicitly incorporates association (kon) and dissociation (koff) kinetics, accounts for rebinding, and extends to systems with drug removal. For drugreceptor closed systems at equilibrium, we prove that relative Ehrlich occupancy time equals b0/(Kd +b0), where Kd is the dissociation constant and b0 the drug concentration. For induced-fit mechanisms, conformational changes reduce the effective dissociation constant to [Formula] (where k3 and k4 are forward and reverse isomerization rates), prolonging occupancy through kinetic trapping. Critically, for drug-receptor systems with first-order drug elimination at rate k3, we derive rigorous bounds: b0/[(b0 + Kd)k3] [≤] EOT{infty} [≤] b0/(Kd {middle dot} k3), revealing that both binding affinity and elimination rate jointly determine occupancy. This explains why high-affinity drugs can fail clinically if eliminated rapidly, and identifies pharmacokinetic optimization opportunities. We prove Copelands definition is a special case of Ehrlich occupancy time when rebinding is absent. Our framework provides quantitative tools for optimizing drug design beyond binding affinity and enables improved prediction of in vivo efficacy where pharmacokinetics dominate.

Since the introduction of vaccination in the current COVID-19 outbreak, many countries have approved and implemented vaccination campaigns to mitigate and ultimately curtail the pandemic. Several types of vaccines have been proposed and many of them have finally been approved and used in different countries. The different types of vaccines have different vaccine parameters, and therefore, this situation induces the necessity of modeling mathematically the scenario of multiple imperfect vaccines. In this paper, we introduce a SIR-based model considering different vaccines, and study the basic properties of the model, including the stability of the Disease-Free Equilibrium (DFE), which is locally asymptotically stable if the reproduction number is less than 1. A sequence of further results aims to enumerate the conditions where the reproduction number can be decreased (or increased). Two important mathematical propositions indicate that in general vaccination might not be enough to contain an outbreak and that the addition of new vaccines could be counterproductive if the leakiness parameter is greater than a threshold{eta} . This model, despite its simplicity, was validated with data of the COVID-19 pandemic in five countries: Israel, Chile, Germany, Lithuania, and Czech Republic, observing that improvements for the vaccine campaigns can be suggested by the developed theory.

We introduce a mathematical model of the cancer-immunity cycle and use it to test several hypotheses regarding the combination, timing, and optimization associated with chemotherapy and immunotherapy dosing schedules in the context of competition and selection pressure. A key conceptual idea is the value of synchronizing the dosing schedules with the fundamental period of the cancer-immunity cycle. The competitors in the population dynamics evolutionary game are the cancer cells, healthy cells, and T-cells, which form a non-transitive rock-paper-scissor chain, mediated by the tumor microenvironment. The chemotherapy and immunotherapy dosing schedules each act as control functions whose timing we synchronize with the fundamental period of the underlying nonlinear dynamical system. With the model, we show among other more detailed results, that chemotherapy and immunotherapy schedules are non-transitive; the best duration of the chemotherapy is around one-quarter of the cancer-immunity cycle, whereas for immunotherapy it is one-half cycle; immunotherapy dosing should preceed chemotherapy dosing. A general conclusion is that optimized timing of the dosing schedules can make up for lower total dose, opening up new possibilities for designing less toxic and more efficacious dosing regimens with drugs currently in use. Obtaining and calibrating more accurate measurements of the cycle-period across patient populations would be an important step in making some of these ideas clinically actionable.

Multiple Myeloma (MM) is a plasma cell cancer that occurs in the bone marrow. A leading treatment for MM is the monoclonal antibody Daratumumab, targeting the CD38 receptor, which is highly overexpressed in myeloma cells. In this work we model drug evasion via loss of CD38 expression, which is a proposed mechanism of resistance to Daratumumab treatment. We develop an ODE model that includes drug evasion via two mechanisms: a direct effect in which CD38 expression is lost without cell death in response to Daratumumab, and an indirect effect in which CD38 expression switches on and off in the cancer cells; myeloma cells that do not express CD38 have lower fitness but are shielded from the drug action. The model also incorporates competition with healthy cells, death of healthy cells due to off-target drug effects, and a Michaelis-Menten type immune response. Using optimal control theory, we study the effect of the drug evasion mechanisms and the off-target drug effect on the optimal treatment regime. We identify a general increase in treatment duration and costs, with varying patterns of response for the different controlling parameters. Several distinct optimal treatment regimes are identified within the parameter space. Author summaryIn this work we investigate a model of Multiple Myeloma, a cancer of the bone marrow, and its treatment with the drug Daratumumab. The model incorporates proposed mechanisms by which the cancer evades Daratumumab by reduced expression of the receptor CD38, which is the drug target and normally abundent in the cancer cells. The model includes an off-target effect, meaning that the drug treatment destroys some healthy cells alongside the targeted cancer cells. Both mechanisms can reasonably be expected to reduce the efficacy of the drug. We investigate the model using optimal control methods, which are used to find the drug dose over time which best balances the financial and health costs of treatment against cancer persistence, according to a specified cost function. We show that this drug resistence and off-target effect prolongs the optimal treatment and increase the burden of both the disease and drug. We analyse the distinct effects of the controlling parameters on each of these costs factors as well as the time course, and identify conditions under which extended treatment is required, with either intermittant treatment or a steady reduced dose. Extended treatment may be indefinite or for a fixed period.

Antimicrobial resistance is a major threat to global health and food security today. Scheduling cycling therapies by targeting phenotypic states associated to specific mutations can help us to eradicate pathogenic variants in chronic infections. In this paper, we introduce a logistic switching model in order to abstract mutation networks of collateral resistance. We found particular conditions for which unstable zero-equilibrium of the logistic maps can be stabilized through a switching signal. That is, persistent populations can be eradicated through tailored switching regimens. Starting from an optimal-control formulation, the switching policies show their potential in the stabilization of the zero-equilibrium for dynamics governed by logistic maps. However, employing such switching strategies, deserve a specific characterization in terms of limit behaviour. Ultimately, we use evolutionary and control algorithms to find either optimal and sub-optimal switching policies. Simulations results show the applicability of Parrondos Paradox to design cycling therapies against drug resistance.

Achieving complex multi-species control objectives is essential for engineering advanced autoregulated biomolecular devices. This paper addresses the problem of robust steady-state tracking for outputs defined as multiplicative combinations of biomolecular species concentrations. We first introduce a control architecture realized via chemical reaction networks that steers the product of two target species concentrations in the controlled network to a prescribed value. A robust stability analysis is provided for closed-loop system families with distinct structural characteristics. The proposed framework is also extended to a more general formulation capable of regulating arbitrary monomial outputs involving multiple species. Numerical simulations of representative examples corroborate the theoretical results and illustrate the effectiveness of our approach.

This work is focused on Ordinary Differential Equations(ODE)-based models of biochemical systems that possess a singular Jacobian manifesting in non-hyperbolic equilibria. We show that there are several classes of systems that exhibit this behavior: a)systems with monomial-type interaction terms and b)systems with linear or nonlinear conservation laws. While models derived from mass-action principles often present with linear conservation laws stemming from the underlying biologic rationale, nonlinear conservation laws are more subtle and harder to detect. Nevertheless, in both situations the corresponding ODE system will contain non-hyperbolic equilibria. While having a potentially more complex dynamics and falling outside of the scope of existing theoretical frameworks, this class of systems can still exhibit adapting behavior associated with certain nodes and inputs. We derive a generalized adaptation condition that extends to singular systems and is compatible with both single-input/single-output and multiple-input/multiple-output settings. The approach explored herein, based on the notion of Moore-Penrose pseudoinverse, is tested on several synthetic systems that are shown to exhibit homeostatic behavior but are not covered by existing methods. These results highlight the role of the network structure and modeling assumptions when understanding system response to input and can be helpful in discovering intrinsic relationships between the nodes.

Oncolytic virotherapy (OVT) represents an innovative and promising therapeutic method for cancer treatment. This approach involves the introduction of oncolytic viruses into the patient, which are engineered to selectively target and lyse tumor cells. Based on previous mathematical modelling of oncolytic viruses, we consider a mathematical model that describes the intricate interactions between the oncolytic virus, cancer cell populations, and the immune system. Our study includes a detailed qualitative and quantitative analysis of the model to explain why, despite their promise, oncolytic viruses alone rarely lead to complete and lasting regression of established tumors in vivo. We use parameter sensitivity analysis to support our findings. Furthermore, we consider the spatial version of the model in the form of a system of reaction-diffusion equations. A travelling wave analysis shows an unexpected phenomenon. The solution components do not necessarily evolve into a single traveling front; rather, they develop into stacked fronts, where each front propagates at a different speed. We give explicit formulas for these different invasion speeds, confirm those through numerical simulations, and discuss their significance for OVT.

This study presents a theoretical and mathematical framework for understanding the dynamical behavior of infectious disease spread using a compartmental modeling approach. The proposed model incorporates memory effects to capture temporal dependencies that are not adequately represented by classical formulations. Qualitative analysis is employed to investigate the stability properties of the system and the role of key mechanisms in shaping long term dynamics. Publicly available surveillance information is used only to illustrate the consistency of the model behavior with observed trends. The results highlight the value of memory based modeling structures for describing complex biological processes and provide a general mathematical perspective for studying epidemic dynamics.

Motivated by curiosities of disease progression seen in the coronavirus pandemic, we analyze a minimalist predator-prey model for the immune system (predator) competing against a pathogen (prey). We find that the mathematical model alone accounts for numerous paradoxical behaviors observed in this and other infections. These include why an exponentially growing pathogen requires an exposure threshold to take hold, how chronic and recurrent infections can arise, and what can allow very sick patients to recover, while healthier patients succumb. We also examine the distinct dynamical roles that specific, "innate," and nonspecific, "adaptive," immunity play, and we describe mathematical effects of infection history on prognosis. Finally, we briefly discuss predictions for some of the effects of timing and strengths of antibiotics or immunomodulatory agents.

This paper introduces the notion of cumulative dose response (cDR). The cDR is the area under the plot of a response variable, an integral taken over a fixed time interval and seen as a function of an input parameter. This work was motivated by the accumulation of cytokines resulting from T cell stimulation, where a non-monotonic cDR has been observed experimentally. However, the notion is of general applicability. A surprising conclusion is that incoherent feedforward loops studied in the systems biology literature, though capable of non-monotonic dose responses, can be mathematically shown to always result in monotonic cDR.

In this work we use mathematical modeling to describe a possible route to the end of COVID-19, which does not feature either vaccination or herd immunity. We call this route self-burnout. We consider a region with (a) no influx of corona cases from the outside, (b) extensive social distancing, though not necessarily a full lockdown, and (c) high testing capacity relative to the actual number of new cases per day. These conditions can make it possible for the region to initiate the endgame phase of epidemic management, wherein the disease is slowly made to burn itself out through a combination of social distancing, sanitization, contact tracing and preventive testing. The dynamics of the case trajectories in this regime are governed by a single-variable first order linear delay differential equation, whose stability criterion can be obtained analytically. Basis this criterion, we conclude that the social mobility restrictions should be such as to ensure that on the average, one person interacts closely (from the transmission viewpoint) with at most one other person over a 4-5 day period. If the endgame can be played out for a long enough time, we claim that the Coronavirus can eventually get completely contained without affecting a significant fraction of the regions population. We present estimates of the duration for which the epidemic is expected to last, finding an interval of approximately 5-15 weeks after the self-burnout phase is initiated. South Korea, Austria, Australia, New Zealand and the states of Goa, Kerala and Odisha in India appear to be well on the way towards containing COVID by this method.

Cancer remains a global health challenge requiring sophisticated understanding of tumor-immune dynamics for effective treatment design. Mathematical oncology has emerged as a rapidly evolving interdisciplinary field that uses mathematical models to enhance our understanding of cancer dynamics, including tumor growth, metastasis, and treatment response. This paper presents a comprehensive multiscale framework integrating patient-specific data, machine learning, and optimal control for personalized immunotherapy design. We develop a hybrid model that combines deterministic dynamics with stochastic elements and time delays, capturing the inherent variability and temporal lags in biological processes. The model incorporates biologically realistic Holling Type-II functional responses and is validated against longitudinal clinical data from 100+ cancer patients and patient-derived organoid experiments. Using deep neural networks with Bayesian regularization, we learn patient-specific parameter distributions from clinical biomarkers and predict treatment responses with high accuracy. Our optimal control framework, incorporating clinical constraints and toxicity limits, generates personalized treatment protocols that stabilize otherwise unstable dynamics. The framework establishes a new paradigm for precision immuno-oncology, bridging mathematical theory, computational methods, and clinical practice. Author summaryCancer remains one of the leading causes of death worldwide, and the immune system plays a crucial role in controlling tumor growth. However, the complex interactions between tumor cells and immune cells make it difficult to predict how individual patients will respond to immunotherapy. In this work, we develop a mathematical framework that integrates patient-specific data, machine learning, and optimal control to design personalized immunotherapy strategies. Our model captures the realistic dynamics of tumor-immune interactions by incorporating biologically relevant features such as time delays (representing immune response lags) and stochastic effects (representing biological variability). Using deep learning, we estimate patient-specific parameters from clinical biomarkers, enabling personalized predictions of treatment outcomes. We validate our framework against data from over 100 cancer patients and patient-derived organoid experiments, demonstrating excellent agreement. Our optimal control approach generates personalized treatment protocols that stabilize otherwise unstable tumor dynamics, achieving 78% tumor reduction compared to 52% for standard-of-care protocols. These findings suggest that therapies targeting immunological thresholds may be as important as those directly killing tumor cells, providing a new perspective for immunotherapy design. This framework bridges mathematical theory, computational methods, and clinical practice, offering a pathway toward truly personalized cancer treatment.

Non-pharmaceutical interventions were very instrumental in the early phases of COVID-19 pandemic. Fortunately, the urgency to control the crisis prompted an accelerated vaccine development process, saving millions of lives globally. Despite these measures, cases of COVID-19 and other variants SARS-CoV-2 are still being reported in different countries and regions. We present a within-host mathematical model of SARS-CoV-2 infection that incorporates target cell dynamics, innate and adaptive immune responses, and vaccine interventions. Analytical results identify the basic reproduction number, R0, as the threshold for infection persistence. Simulations show that immune responses, both lytic and non-lytic, are critical in controlling viral replication, with vaccine-induced immunity further reducing viral load and protecting epithelial cells. Immune-boosting strategies and monoclonal antibody therapies targeting intracellular replication outperform entry-blocking interventions alone, while combination approaches yield the greatest reduction in peak viral load and fastest clearance. Timing is crucial: early vaccination or treatment maximizes benefits, whereas delays allow higher viral titers to persist. These results underscore the importance of early, multi-mechanism interventions, particularly for vulnerable populations such as the elderly and immunocompromised. The model offers a framework for evaluating treatment strategies and can be extended to incorporate pharmacokinetics/pharmacodynamics or patient-specific calibration for improved predictive accuracy.

Auxins are a family of plant hormones involved in various processes across plant tissues and species. The Nuclear Auxin Pathway (NAP) consists of interacting transcription factors (ARFs) and repressors (Aux/IAAs), which govern an individual cells response to changes in auxin concentration. These components are present in all land plants, and many species possess multiple copies of each signalling component. We present a general framework for ODE-based models of NAP submodules with the flexibility to model the promotion and repression of target genes by any combination of transcriptional regulators. We analyse published data and show that auxin treatment in Arabidopsis thaliana roots triggers a range of characteristically distinct temporal response profiles--for both target genes and the signalling components themselves. Using our modelling framework, we recapitulate aspects of this behaviour by presenting examples of real and theoretical NAP subnetworks, and by analysing the effect that these network dynamics have on auxin-mediated transcriptional responses. This work demonstrates the utility of our modelling framework as a general-purpose tool for understanding the function of certain protein-protein and protein-DNA interactions through their effects on the NAP. This exploration of the rich dynamics of more complex signalling pathways promises to advance our understanding of the NAP.

Type 1 diabetes (T1D) is an autoimmune disease in which the immune system attacks pancreatic beta cells, leading to high blood glucose levels and requiring lifelong insulin therapy. There is no cure, and individuals with T1D may face a reduced lifespan of up to 12 years. Defects in regulatory T cells (Tregs) are a key contributor to disease onset and are being explored as a therapeutic avenue. However, the effectiveness of Treg therapy remains uncertain. Research is further limited by the inability to directly observe pancreatic and lymph node activity during the long presymptomatic stage of T1D. In this study, we develop a mathematical model for beta and T cell dynamics. We find both Treg quality and quantity affect disease progression, and that antigen-presenting cell (APC) dynamics play a central role. Notably, Treg therapy combined with APC depletion improves outcomes, especially with strong peptide-induced APC activation.